#include <stdio.h>
#include <time.h>
#include <limits.h>
#include <gmp.h> // gcc -lgmp

#define N 25
#define BASE 10
#define MANTISSA 1024*1024*32

/*
N = 10.000
MANTISSA de 128K
104.00 segundos

*/

void a_zero(mpf_t rop){ // OK = 0.34314...

	mpf_t aux;
	mpf_init(aux);
	
	// a = 6.0-(4.0*sqrt(2));
	mpf_sqrt_ui(aux, 2); 
	mpf_mul_ui(aux, aux, 4);
	mpf_ui_sub(rop, 6, aux);
	
	mpf_clear(aux);
}

void y_zero(mpf_t rop){ // OK = 0.41421...

	mpf_t aux;
	mpf_init(aux);
	
	//y = sqrt(2)-1.0;	
	mpf_sqrt_ui(aux, 2); 
	mpf_sub_ui(rop, aux, 1);

	mpf_clear(aux);
}

void y_kn(mpf_t rop, mpf_t y){ // Ok 0.3734i...^e-2

	mpf_t y4, n, d;	
	mpf_init(y4);	
	mpf_init(n);
	mpf_init(d);
	
	//	yn = (1-(1-y^4)^1/4)/
	//	     (1+(1-y^4)^1/4)

	mpf_pow_ui(y4, y, 4);
	mpf_ui_sub(y4, 1, y4);
	mpf_sqrt(y4, y4);
	mpf_sqrt(y4, y4);
	
	mpf_ui_sub(n, 1, y4);

	mpf_add_ui(d, y4, 1);

	mpf_div(rop, n, d);

	mpf_clear(y4);
	mpf_clear(n);
	mpf_clear(d);
}

void a_kn(mpf_t rop, mpf_t a, mpf_t yn, unsigned long k){
	
	mpf_t a4, exp, yk;

	//an = a(1+ yn)^4 - 2^(2k+3) *  yn(1 + yn + yn^2)
	//**	an converge para 1/pi
	
	// a4 = a(1+yn)^4
	mpf_init(a4);
	mpf_add_ui(a4, yn, 1);
	mpf_pow_ui(a4, a4, 4);
	mpf_mul(a4, a, a4);

	// exp = 2^(2k+3)
	//mpf_init2(exp, MANTISSA*2); // Pode virar um numero muito grande
	mpf_init(exp);
	mpf_set_ui(exp, 2);
	k *= 2;
	k += 3; 
	mpf_pow_ui(exp, exp, k);	

	// yk = yn(1+yn+yn^2)
	mpf_init(yk);
	mpf_pow_ui(yk, yn, 2);
	mpf_add(yk, yk, yn);
	mpf_add_ui(yk, yk, 1);
	mpf_mul(yk, yk, yn);	

	// an = a4-exp*yk
	mpf_mul(exp, exp, yk);
	mpf_sub(rop, a4, exp);

	mpf_clear(a4);
	mpf_clear(yk);
	mpf_clear(exp);
}

int main(int argc, char* argv[]){

	unsigned long int k;
	time_t t_zero, t_final;
	double delta_t;

	time(&t_zero);

	mpf_set_default_prec(MANTISSA);
	mpf_t a, y, pi;

	mpf_init(a);
	a_zero(a);

	mpf_init(y);
	y_zero(y);
	
	/*

	a = 6 - (4 * 2^1/2)
	y = (2^1/2) - 1

	yn = (1-(1-y^4)^1/4)/
	     (1+(1-y^4)^1/4)

	an = a(1+ yn)^4 - 2^(2k+3) ? yn(1 + yn + yn^2)
	**	an converge para 1/pi

	*/
	
	for(k=0; k<N; k++){// K na primeira interação é = 0
		y_kn(y, y);
		a_kn(a, a, y, k);
    }

	mpf_clear(y);

	mpf_ui_div(a, 1, a); // an converge para 1/pi
	mpf_out_str(NULL, BASE, 0, a);
	mpf_clear(a);
	
	time(&t_final);
	delta_t = difftime(t_final, t_zero);
	printf("\nTempo = %.2lf segundos\n", delta_t);	

	return 0;
}
